Climatic control on the location of continental volcanic arcs

Orogens and volcanic arcs at continental plate margins are primary surface expressions of convergent plate tectonics. Although it is established that climate affects the shape, size, and architecture of orogens via orographic erosion gradients, the ascent of magma through the crust and location of volcanoes along magmatic arcs have been considered insensitive to erosion. However, available data reveal westward migration of late-Cenozoic volcanic activity in the Southern Andes and Cascade Range where orography drives an eastward migration of the topographic water divide by increased precipitation and erosion along west-facing slopes. Thermomechanical numerical modeling shows that orographic erosion and the associated leeward topographic migration may entail asymmetric crustal structures that drive the magma ascent toward the region of enhanced erosion. Despite the different tectonic histories of the Southern Andes and the Cascade Range, orographic erosion is a shared causal mechanism that can explain the late-Cenozoic westward migration of the volcanic front along both magmatic arcs.

www.nature.com/scientificreports/ posed where break-up between the JDF and the Explorer plates occurs accommodated by the transform Nootka Fault entering into subduction 21,27 , and along the southern edge of the JDF plate, nearby the Mendocino Triple Junction 19,21 . Emplacement of volcanic and plutonic rocks in the North Cascades (48-50°N) occurs through compressive and transcurrent shear zones, such as the Straight Creek and Ross Lake fault, that accommodate the accretion of allochthonous terranes since the Mesozoic [28][29][30] . These shear zones were suggested as preferential conduits for the emplacement of the Chilliwack Batholith, which crystallized from ~ 35 to ~ 4 Ma and is considered of the same magmatic origin as the Quaternary Mt. Baker volcano 23,31-34 (transect Fig. 1a, b). The volcanic arc in the North Cascades migrates from the northeast towards the southwest since at least ~ 4 Ma. If one considers the intrusive rocks of the Chilliwack Batholith the onset of such trend can be dated back to ~ 12 Ma (Refs. 23,[31][32][33]. Conversely, in the High Cascades (40-48°N) volcanic rocks erupted since the late Miocene are preserved, providing the record of an eastward migration of the arc front 23,[31][32][33][34] in a predominantly transtensional tectonic regime that extends into the Basin and Range sector 19,20,34,35 .
The North Cascades interact with westerly Pacific winds, which generate ten times higher precipitation rates on western slopes with respect to the eastern side of the orogen at latitudes where the westerlies are perpendicular to the orogenic belt 36 (north of ~ 46° N, Fig. 1a). This precipitation pattern is consistent with focused rock exhumation in the western flank of the North Cascades since ~ 10 Ma (Fig. 1b), suggesting that orographic effects are long term features in western North America 11,14,[37][38][39][40] . While Mesozoic to Tertiary plutonic rocks and metamorphic basement are exhumed in the North Cascades due to orographic erosion, Tertiary to Quaternary volcanics are still preserved in the High Cascades 23,30,31,33,37 . Accordingly, the topographic water divide north of ~ 46°N is 40-70 km more to the east and the average elevation is lower by about 500-1000 m than in the southern sector ( Fig. 1a, b). The water divide migration is supposed to have occurred during the Plio-Quaternary 14 . Quaternary volcanoes north of ~ 46°N commonly overlie slab isodepths < 100 km and are systematically located west of the topographic water divide (Fig. 1a). Quaternary volcanoes south of ~ 46°N, instead, overlie slab isodepths ≥ 100 km and are found along the topographic water divide.
The Southern Andes Volcanic Zone (SVZ) is generated by the subduction of the Nazca oceanic plate below the western margin of the South American continent between latitudes 33-46°S (Fig. 1c). Subduction occurs at an average rate of ~ 7 cm/yr (Ref. 41 ) with an average slab dip angle of ~ 25° (Refs. 42,43 ) since the late Miocene 44 . The age of the subducting oceanic plate decreases southward from ~ 40 Ma to present-day where the Chile Ridge is currently entering into subduction at the Chile Triple Junction (CTJ, Ref. 45 ). In the southern Central Andes (33-38°S) several fold-and-thrust belts accommodate deformation and volcanic emplacement since at least the Miocene 44,46,47 . In the northern Patagonian Andes (38-46°S) the general compressive tectonic regime changes to transpressive accommodated by the major strike-slip Liquiñe-Ofqui Fault Zone (LOFZ) active since at least ~ 5 Ma, with denudation ages between ~ 16-10 Ma suggesting an earlier activation [48][49][50][51] . In this sector, the Patagonian Batholith is exposed, which allows for reconstructions of the magmatic arc migration since the Cretaceous 52 . Conversely, batholiths are not exposed in the Central Andes where deformed volcanic and sedimentary rocks are preserved at the surface 44,46 . Eruption and emplacement ages of volcanic and intrusive rocks south of 40°S suggest a westward migration and narrowing of the southern volcanic front by around 50 km into the region of the LOFZ since at least the Pleistocene 47,53-55 (e.g., Tronador-Osorno transect, Fig. 1c, d).
Similarly to the Cascade Range, south of 40°S the SVZ interacts with westerly Pacific winds perpendicular to the orogenic belt associated with significantly increased precipitation on the western slopes of the Southern Andes 36 (Fig. 1c). Exhumation ages younger than ~ 10 Ma are recorded in the western part of the orogen, whereas systematically older exhumation ages are recorded on the drier eastern peaks of the orogen 49,50,56 (Fig. 1d). An abrupt increase in exhumation rates since ~ 7 Ma is recorded in the Patagonian Andes, pinpointing the onset of Andean glaciations 40,49,[56][57][58] . An erosion hotspot at 42-46°S is observed in thermochronologic data 2 ± 2 Ma (Fig. 1c, Ref. 58 ). In addition, an overall decrease in mean elevation from about 3 km to 1 km and an ~ 70 km eastward migration of the topographic water divide further testify the long-term regional orographic effects on the Andean topography 56,59 . Quaternary stratovolcanoes overlie slab isodepths ≤ 100 km south of 40°S and are systematically west of the topographic water divide (Fig. 1c) whereas, north of 40°S, they overlie slab isodepths > 100 km and are found along the main topographic divide.
The observed topographic leeward and volcanic arc front windward migrations, shared by both the CVA and SVZ despite the different tectonic history, suggest a potential climatic control via erosion on the long-term evolution and location of these continental volcanic arcs. Hereafter, we assess the plausibility of a coupling between crustal magma transfer and climate-controlled erosion.
Coupled geodynamic and erosion numerical modeling. We use a thermo-mechanical (visco-elastoplastic) geodynamic numerical model coupled to the stream power erosion model (see Methods) to test the hypothesis of an upwind (i.e., westward) migration of the magma ascent and volcanic arc front forced by orography-driven topographic change. The model accounts for the continental crust, upper mantle, asthenosphere (rheological properties are shown in Supplementary Table 2) and a wedge-shaped topography up to 2.5 km high. In a first set of numerical experiments, we impose the initial topography directly above or shifted to the right with respect to a central mantle melting zone (MMZ), analogue to the region where partial melting is generated at ~ 100 km depth feeding crustal magma reservoirs and surface volcanoes 1-4 ( Fig. 2 and Supplementary  Figs. 1-4), and no erosion is imposed. This first set of numerical experiments allows us to assess changes in the magma ascent prior to and after the topographic landscape has adjusted to an orographic erosional gradient. In a second set of simulations we impose asymmetric surface erosion, with stream power erosion rates being one order of magnitude higher on left-facing slopes than on right-facing slopes. In these experiments, we consider an initial topography directly above the MMZ ( Fig. 3  www.nature.com/scientificreports/  www.nature.com/scientificreports/ continuously adjusting to the progressive crustal deformation due to the emplacement of magma, modulate the erosion rates 11,61 (see Eq. 13 in Methods), allowing us to address the transient evolution of magma ascent as the topography adjusts to orographic erosion. In order to isolate the effects of climate-driven orography over the crustal trajectories of magma emplacement and ascent towards the surface, we did not impose any tectonic convergence, subduction of oceanic plate, crustal accretion nor sedimentation in these experiments, which simply evolve by upwelling of buoyant magma into the crust until depletion of the MMZ over about 0.3 Ma. This is obviously an oversimplification, but it reduces the degree of complexity of the model and removes any effect from poorly constrained parameters (e.g. spatial-temporal changes in convergence rate or thermal/mechanical/ rheological properties). Although this approach does not allow to reproduce the geologic history of the CVA and SVZ in our simulations, it enables us to explore plausible orographic effects on the joint evolution of the topographic landscapes, crustal structures and magma ascent observed in these settings, that is our ultimate goal. We obviously bear this limitation in mind when framing our results in the contexts of the CVA and SVZ. The model domain extends over 200 × 120 km (horizontal, x, and vertical, y, dimensions), resolved by 201 × 61 grid points respectively, and is distributed on an irregular Eulerian grid that accounts for a maximum resolution of 1 km along the x and y directions in the upper-central part of the model domain. 400 × 400 Lagrangian markers are randomly distributed along the x and y dimensions and used for advecting the material properties 16,62,63 (Supplementary Table 2). The material properties carried by Lagrangian markers are then interpolated onto the Eulerian grid via a 4th order Runge-Kutta interpolation scheme 62,63 . An internal free surface is simulated through a 10 -km thick layer of sticky air 62 . The velocity boundary conditions are free slip at all boundaries (x = 0 and x = 200 km; y = 0 and y = 120 km). The initial temperature gradient is piece-wise linear resulting from an adiabatic temperature gradient of 0.5 °C/km in the asthenosphere 64 and thermal boundary conditions fixed at 0 °C at the surface and 1327 °C at the lower boundary, with nil horizontal heat flux across the vertical boundaries 63 . Temperature is set to 1327 °C in the MMZ, that is initially circular with 20 km diameter and imposed in the center of the model domain at 100 km depth, consistent with estimates from the CVA 21,65 and SV 1,66,67 . The upward transfer and emplacement of magma into the crust occurs by buoyancy through a 3 -km wide magmatic channel, which allows for a simplified representation of magmatic percolation by hydrofractures, diffusion, porous flow, and reactive flow through the rheologically stronger mantle lithosphere 63 . The bulk composition of the MMZ www.nature.com/scientificreports/ is ultramafic, but enrichment by 25% mafic melt is prescribed in the magmatic channel 63 . The viscosity, η ductile (Eq. 3 in Methods), of partially molten rocks (with ξ > 0.1, Eq. 10 in Methods) is assigned a low constant value of 10 16 Pa s, and the η ductile upper limit is set to 10 25 Pa s (Refs. 16,63,64 ). The topographic wedge measures 2.5 km in simulations with central topography, and 1.5 km in simulations with lateral topography, consistent with the observed orography-related change in mean altitude in the CVA and SVZ (Refs. 38,39,56,59 ). The parametric study focuses on three main parameters. First, the crustal thickness (Ch) varies between 45 and 35 km, consistently with literature values 42,65 . Second, the initial depth of the magmatic channel upper tip (Dmc), which assumes values of 15, 20, and 25 km to account for efficient/inefficient magma transfer through the crust. Values are based on estimates of the depth of the regional brittle-ductile transition (~ 20 km, Refs. 42,51,65 ), assuming that the rheology of the upper plate crust exerts a primary control on the efficiency of magma transfer. Third, the thickness of the thermal lithosphere (Lh) varies between 90 and 100 km allowing us to account for variations of the elastic lithosphere between ~ 55 and ~ 65 km, as estimated from the regional geophysical data 42,65 .

Results
Results from the first set of numerical experiments show that the emplacement of magma into the crust is accommodated by structures that affect both the brittle upper crust and the ductile lower crust. Simulations that account for an initial central topography allow the formation of symmetric crustal magma reservoirs and structures with respect to the MMZ, feeding volcanoes on both sides of the topographic wedge (Fig. 2a, b and Supplementary Figs. 1-4). Simulations that account for an initial lateral topography develop highly asymmetric crustal magma reservoirs and structures, feeding volcanoes on the opposite side of the initial topographic wedge (Fig. 2c, d Table 1 and Fig. 4). For all tested configurations, a thin crust (Ch = 35 km) and a cold lithosphere (Lh = 100 km) facilitate the brittle strain ( Fig. 2 and Supplementary Figs. 1, 2). When the initial topography is lateral, the leftward asymmetry of the strain and magma ascent is enhanced (Figs. 2c-d, 3, and Supplementary Figs. 1, 2). When the initial depth of the magmatic channel upper tip is deep (Dmc = 25 km), the viscous strain in the lower crust accommodates the emplacement of magma, inhibiting the formation of brittle upper crustal structures and the rise of magma to the near surface ( Supplementary Figs. 1-4). When Dmc is shallow (15 km), the magma rises easily to the surface through brittle structures near the center of the model, resulting in relatively small values of δd ( Fig. 4 and Supplementary Figs. 1-4). When Dmc is intermediate (20 km), brittle-ductile structures cross through the crust reaching the surface at further distance from the center of the model (Figs. 2-3, and Supplementary Figs. 1, 3, 4). In the reference model (i.e., Ch = 35 km, Lh = 100 km, and Dmc = 20 km, Fig. 2), an initial lateral topography leads to δd ≈ −45 km (i.e., leftward asymmetry), comparable to the observed westward migration of the magmatic arcs in the natural case studies (Fig. 1). We further remark that, in simulations with initial central topography, a symmetric set of thrusts rooted in the region of magma emplacement form a pop-up structure similar to those observed in the northern sector of the SVZ 46,47 . In the simulations with initial lateral topography, the asymmetric thrust verging towards the west is similar to the "Western Patagonian Thrust" proposed by Ref. 51 .
Results Uplift of additional topographic peaks on the left side of the model domain occurs due to left-verging thrusts accommodating the magma ascent. We associate this topography to that of volcanic edifices west of the main topographic divide (Fig. 1). When the crust is thick (Ch = 45 km) the viscous strain accommodates almost entirely the magma emplacement and limited bulk strain accumulates along shallow brittle faults that tend to be symmetric ( Fig. 3 and Supplementary Figs. 7, 8). Accordingly, no topographic relief is generated on the left side of the model domain and the topography migrates to the right by ~ 10 km in ~ 0.3 Ma (Supplementary Figs. 7, 8).

Discussion
The break-up of the Farallon plate into the Juan de Fuca and Nazca plates in the early Miocene established a kinematic configuration similar to the present day in terms of slab partitioning and convergence rates along the American Cordillera 20,44 . Although the onset of the eastward migration of the topographic water divide in the CVA and the SVZ is uncertain, systematically younger exhumation ages on the western slopes of both the North Cascades and northern Patagonian Andes (Fig. 1) suggests that orographic effects were already established in the latest Miocene 37-39,56,57 , leading to the formation of the current landform throughout the Pleistocene 14,39,58,59 . Although low-temperature thermochronometric ages constrain the onset of enhanced erosion between ~ 10 and ~ 6 Ma in the North Cascades and the Southern Andes, respectively 37,38,56 , large uncertainties exist about the paleo position of topographic water divide and paleo-elevations in general. A similar issue pertains to the onset of arc front migration, which can be only roughly constrained from batholith ages 31,52 , but precise information about the exact position of the paleo volcanic arc front is lacking. Our modeling results show a fast response of the crustal strain and magmatic plumbing system to orographic erosion, which adjust to topographic changes in just a few hundreds of thousands of years (Figs. 2 and 3). www.nature.com/scientificreports/ In the CVA, the present-day clockwise rotation of the Oregon and Washington blocks and the counterclockwise displacement of the volcanic front since the Miocene have been ascribed to differential along strike rollback of the Juan de Fuca slab 20,31,34,35 . However, since the magmatic source follows the rollbacking slab 1,4-6 , this model requires some unjustified degree of decoupling between rollback and arc front migration 34 . In the North Cascades, the increase in the rate of arc front migration towards the southwest since ~ 4 Ma (Refs. [31][32][33] ) requires faster slab rollback rates than during the late Miocene, difficult to reconcile with the transpressive regime and present northeastward displacements of the Washington block 31,34 . The presence of slab-derived magmas (i.e., adakites 68 ) in southernmost and northernmost volcanoes of the CVZ may suggest shallower (i.e., closer to the trench) magma production near the breakup of the Juan de Fuca Plate into the Gorda and Explorer plates 22,24,68 . However, if the thermal gradient along the strike of the subduction system set alone the surface location of the volcanic arc, trench-ward migration of the volcanic front would be observed both in the northern and southern sectors of the CVA. Thus, the opposite sense of arc front migration in the North Cascades and in the High Cascades 34 can hardly be explained by along-strike slab age and/or thermal changes alone. Since Quaternary volcanoes are emplaced over fault systems, orographic interactions with westerly Pacific winds (Fig. 1a, b) and the mechanistic link between asymmetric erosion and crustal structures accommodating the magma ascent (Figs. 2, 3, 4 and 5) appear as a suitable along-strike differential forcing to drive the observed latitudinal arc migration front trend reversal. Magnetotelluric data (Fig. 5) show westward dipping magmatic plumbing systems in the North Cascades 69 , consistent with our modeling results (Figs. 2c, d and 3) and an orographic control on the location of volcanic arc front.
In the Southern Andes, a Plio-Quaternary magmatic arc narrowing toward the west has been related to the steepening of a formerly shallow slab 55 or a decrease in convergence velocity 54 . South of ~ 40°S, Plio-Pleistocene strain localization along the transpressive system of the LOFZ has been ascribed to oblique convergence of the Nazca Plate and the CTJ approaching its current position [48][49][50] , driving the ascent of magmas and setting the observed trench-ward deflection of the SVZ 8,48,49,54 (Fig. 1c). Geophysical imaging of the deep crust and mantle show that the crustal magma reservoirs of the Osorno volcano are located up to 10 km east of the volcanic edifice 70 , consistent with our modeling results and hypothesis of focused orographic erosion facilitating the strain and magma ascent across the windward side of orogens [11][12][13] (Fig. 5). Field evidence regarding the presence of the west vergent "Western Patagonian Thrust" 51 reinforces our proposal of asymmetric faults assisting magma ascent towards the west. Also in agreement with previous studies 8 , our models indicate that orographic erosion of a thin crust facilitates the formation of major structures, suggesting that orographic erosion may have contributed to the activation, or re-activation, of the LOFZ. Heat advection from the younger slab and the CTJ in the south can likely enhance this effect. The Nazca Plate is younger toward its southern edge (CTJ) and thus the thermal regime of the subduction is likely higher in the northern Patagonian Andes 53,67 . However, no slabderived magmas are found in the whole SVZ 8,52,66,71 and the depth of the source of partial melting is unlikely to change significantly along the strike of the belt 1,52,67 . Differences in magma composition along the strike of the SVZ show a tendency toward more mafic magmas in the south. This trend can be simply explained by a reduced crust-magma interaction due to the thinner crust in this sector 8,66,71 , which also promotes the orographic forcing on the magma ascent proposed here (Fig. 4).  35  90  15  Central  No  65  60  -5   35  90  20  Central  No  70  65  -5   35  90  25  Central  No  10  10  0   35  90  15  Lateral  No  50  10  -40   35  90  20  Lateral  No  50  5  -45   35  90  25  Cateral  No  64  64  0   35  100  15  Central  No  15  15  0   35  100  20  Central  No  65  65  0   35  100  25  Central  No  65  65  0   35  100  15  Lateral  No  15  5  -10   35  100  20  Lateral  No  45  0  -45   35  100  25  Lateral  No  50  We propose a common explanation to the westward migration of the volcanic arc fronts in the North Cascades and the northern Patagonian Andes, located several thousands of kilometers apart and in different tectonic/ structural settings, by the shared orographic interactions between topography and westerly Pacific winds (Fig. 5). The climatic control on magmatism would occur not only through a forcing from ice building-melting and erosion on the magma production 15,72 , but also through orography and topographic change facilitating crustal strain and magma transfer verging toward the region of enhanced precipitation and erosion. Because volcanic arcs provide a substantial contribution to the evolution of climate across timescales, this recognition provides additional evidence of the tight coupling between climate, surface processes, magmatism, and plate tectonics 73

Mechanical component. The continuity equation,
where ρ is the local density, t is time, v is the velocity vector, and ∇ is the divergence operator, allows for the conservation of mass during the displacement of a geological continuum.
(1) ∂ρ ∂t + ∇(ρv) = 0, where η diff and η disl are the shear viscosity for diffusion and dislocation creep, respectively, η 0 is the material static viscosity, σ cr is the diffusion-dislocation transition critical stress, n is the stress exponent, E a is the activation energy, V a is the activation volume, P is pressure, R is the gas constant, T is temperature, and ε II is the second invariant of the strain rate tensor. Then, the viscous deviatoric strain rate tensor, ε ′ ij(viscous) , is computed as: where σ ′ ij is the deviatoric stress tensor,δ ij is the Kronecker delta, ε kk is the volumetric strain rate (e.g., related to phase transformations), η bulk is the bulk viscosity.
Elastic deformation is reversible and assumes proportionality of stress and strain. The elastic deviatoric strain rate tensor, ε ′ ij(elastic) , is computed as: where µ is the shear modulus and ⌣ Dσ ′ ij Dt is the objective co-rotational time derivative of the deviatoric stress tensor. Plastic (or brittle) localised deformation occurs at low temperature in the upper part of the lithosphere, after reaching the absolute shear stress limit, σ yield , is defined as: where C is cohesion and ϕ is the effective internal friction angle. The plastic strain rate tensor, ε ′ ij(plastic) , is then computed as: where X is the plastic multiplier which satisfies the plastic yielding condition σ II = σ yield .
At the lithospheric scale, all deformation mechanisms occur jointly and the overall visco-elasto-plastic rock strain rate tensor, ε ′ ij(bulk) , is defined as: σ yield = C + sin (ϕ)P, www.nature.com/scientificreports/ Partial melting component. Partial melting occurs between the wet solidus, T s , and dry liquidus, T l , of the considered lithologies [75][76][77] (Supplementary Table 2). ξ is the volumetric fraction of melt that increases linearly with T at a constant P (Ref. 64 ), so that, The effective density, ρ eff , of partially molten rocks is then computed by 16,63 : where ρ 0 l and ρ 0 s are the standard densities of solid and molten rocks, respectively. The solid rock density, ρ s , is calculated as 63 : where β is compressibility, α is thermal expansion and ρ o , P 0 , and T 0 are density, pressure and temperature of rocks at surface conditions. Erosion model. During the fluvial incision of an uplifting topography the erosion rate, ė , depends primarily on channel slope, river discharge and rock type 61 . An empirical relationship was derived 11,61 , such that: where z is the elevation, A is the basin drainage area taken as a proxy for river discharge, dz dx is the magnitude of the local topographic slope, and k , m , and n are empirical parameters, usually determined by fitting models to river longitudinal profiles. While k is used as a scaling constant, different values of m and n lead to differences in the shape of river profiles, but commonly used values (e.g., 0.3 < m < 0.5 and n ≈ 1) predict realistic concave-up profiles such as those observed along bedrock rivers incising uplifting topographies 11,61,78 . Although simplistic, Eq. (13) captures the main physical processes and parameter dependencies and, most important for this study, it provides the feedback mechanisms between tectonic uplift and erosion. Tectonic uplift increases slope, and the rise of topography increases drainage area, both leading to an increase in erosion rates.
Integration on a discrete topography of the top of the lithosphere allows for Eq. (13) to be solved numerically, thereby assessing surface elevation changes in response to the tectonic strain and fluvial incision 11,16 . At each time step, the surface load changes associated with modifications of the modelled landscape are computed. The orographic enhancement of precipitation, and hence fluvial erosion, is included in the surface processes model by assigning different k values (Eq. 13) on negatively (left-facing) and positively (right-facing) sloping topography.

Data availability
Supplementary Tables and Figures are available for this paper.